Mapas con R

Visualización de datos geoespaciales con R.

En este taller veremos como visualizar datos geoespaciales (vectoriales) con R, sin embargo, los paquetes de R permiten hacer mucho más que sólo visualizar. Para aprender más sobre esto ver Geocomputation with R (Lovelace, Nowosad, y Muenchow 2019).

💻Preparación

Instala los paquetes y descarga el material del taller.

install.packages(c("usethis", "tidyverse", "sf", "leaflet"))
usethis::use_course("https://github.com/tereom/taller-mapas-mv/blob/master/material.zip?raw=true")

🗺Estáticos

Usaremos los paquetes tidyverse (Wickham et al. 2019) y sf(Pebesma 2018).

library(tidyverse) # arreglar datos (dplyr) y graficar (ggplot2)
library(sf) # simple features para leer y manipular datos espaciales

🐶 Nuestro primer mapa

• Objetivo: Un mapa de población por departamento de Uruguay

• ¿Qué necesito?

  1. Datos que determinan la forma de los departamentos (shapefile),
     
  2. Datos de población por departamento (csv)

Los datos geoespaciales de los departamentos los obtenemos de mapas vectoriales del INE y usaremos la función st_read() para leerlos:

depart_uy <- st_read("datos/ine_depto")

Reading layer `ine_depto' from data source `/Users/tereom/Documents/mapas-taller-mv/datos/ine_depto' using driver `ESRI Shapefile'
Simple feature collection with 20 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
Projected CRS: WGS 84 / UTM zone 21S

La función st_read() recibe como argumento la ruta a la carpeta (puede ser al archivo con terminación .shp) y la asigna a depart_uy.

Los podemos graficar directo con plot(depart_uy) o con ggplot (❤️) en el segundo caso usamos geom_sf() para agregar la capa de coordenadas:

ggplot(depart_uy) +
  geom_sf()

Si inspeccionamos el objeto depart_uy descubrimos que es un poco familiar (y un poco nuevo), lo familiar es que es una tabla de datos de tipo data.frame, con una columna adicional de geometría que almacena la información de la geometría espacial.

class(depart_uy)

[1] "sf"         "data.frame"

glimpse(depart_uy)

Rows: 20
Columns: 6
$ AREA_KM2_  <int> 530, 11928, 4536, 6106, 11643, 10417, 10016, 1392…
$ PERIMETER  <dbl> 202136.96, 722385.94, 404534.80, 431056.43, 82132…
$ DEPTO      <int> 1, 2, 3, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 17, …
$ NOMBRE     <chr> "MONTEVIDEO", "ARTIGAS", "CANELONES", "COLONIA", …
$ CDEPTO_ISO <chr> "UYMO", "UYAR", "UYCA", "UYCO", "UYDU", "UYFD", "…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((579197 6136..., MULT…

En la gráfica podemos colorear (como en cuaqluier ggplot) los polígonos utilizando con el valor de alguna de las columnas del data.frame.

ggplot(depart_uy) +
  geom_sf(aes(fill = AREA_KM2_))

También podemos unirle una tabla de datos externa, por ejemplo de población (datos de juntas departamentales) como lo haríamos de normal:

basicos <- read_csv("datos/datosbasicosjds.csv")
glimpse(basicos %>% select(1:10, 38))

Rows: 19
Columns: 11
$ `Junta Departamental` <chr> "Artigas", "Canelones", "Cerro Largo",…
$ `Tiene WEB`           <chr> "NO", "Si", "Si", "Si", "Si", "Si", "S…
$ `Link de la web`      <chr> NA, "http://www.juntadecanelones.gub.u…
$ Dirección             <chr> "Cnel. Carlos Lecueder 510", "Luis Alb…
$ `Teléfono 1`          <chr> "4772 4833", "4332 2420", "4642 2283",…
$ `Telefono 2`          <chr> NA, NA, "4642 3471", NA, "4362 4778", …
$ Email                 <chr> "juntadeartigas@gmail.com", "juntacan@…
$ Plenarios             <chr> "1er - 2do y 3er jueves de cada mes", …
$ `Hora Plenario`       <chr> "20.30", "18.00", "19.00", "20.00", "1…
$ `Local propio`        <chr> "Si", "Si", "Si", "Si", "Si", "No", "S…
$ Población             <chr> "73.377", "520.173", "84.698", "123203…

# un poco de limpieza para poder unir las tablas
depart_pob <- basicos %>% 
  select(nombre = `Junta Departamental`, pob = `Población`) %>% 
  mutate(nombre = janitor::make_clean_names(nombre), 
         pob = as.numeric(str_remove_all(pob, "\\.")))

depart_uy_pob <- depart_uy %>% 
  mutate(nombre = janitor::make_clean_names(NOMBRE)) %>% 
  left_join(depart_pob)

glimpse(depart_uy_pob)

Rows: 20
Columns: 8
$ AREA_KM2_  <int> 530, 11928, 4536, 6106, 11643, 10417, 10016, 1392…
$ PERIMETER  <dbl> 202136.96, 722385.94, 404534.80, 431056.43, 82132…
$ DEPTO      <int> 1, 2, 3, 5, 6, 8, 9, 11, 12, 13, 14, 15, 16, 17, …
$ NOMBRE     <chr> "MONTEVIDEO", "ARTIGAS", "CANELONES", "COLONIA", …
$ CDEPTO_ISO <chr> "UYMO", "UYAR", "UYCA", "UYCO", "UYDU", "UYFD", "…
$ nombre     <chr> "montevideo", "artigas", "canelones", "colonia", …
$ pob        <dbl> 1318755, 73377, 520173, 123203, 57084, 67047, 588…
$ geometry   <MULTIPOLYGON [m]> MULTIPOLYGON (((579197 6136..., MULT…

Utilizamos el paquete janitor para arreglar los nombres y que coincidan entre las dos tablas que vamos a unir. Unimos con la función left_join() como haríamos en una manipulación de datos ususal de data.frames.

Y graficamos nuestro mapa coloreado por población.

ggplot(depart_uy_pob) +
  geom_sf(aes(fill = pob))

😕 Lo no tan familiar

Shapefiles

Ahora inspeccionamos el lado nuevo de nuestro objeto, cuando leímos los datos apuntamos a una carpeta:

depart_uy <- st_read("datos/ine_depto")

esta carpeta de donde leímos los datos tiene muchos archivos.

fs::dir_tree("datos/ine_depto/")

datos/ine_depto/
├── ine_depto.dbf
├── ine_depto.prj
├── ine_depto.sbn
├── ine_depto.sbx
├── ine_depto.shp
├── ine_depto.shx
└── ine_depto_metadatos.pdf

Esto es porque, en este ejemplo, nuestros datos vectoriales están almacenados en un shapefile (un formato muy común para datos geoespaciales).

Un shapefile es un grupo de archivos que contienen geometrías e información de acuerdo a un estándar especificado por el Insituto de Investigación en Sistemas de Ecosistemas (ESRI). En nuestro ejemplo tenemos los siguientes:

• nombre.shp: es el archivo principal y contiene la geometría.

• nombre.dbf: es la base de datos asoiciada al objeto y almacena la información de los atributos de los objetos, en el caso de ine_depto hay un renglón por polígono, es decir, cada línea contiene las características de un departamento.

ine_depto_dbf <- foreign::read.dbf("datos/ine_depto/ine_depto.dbf")
head(ine_depto_dbf)

  AREA_KM2_ PERIMETER DEPTO     NOMBRE CDEPTO_ISO
1       530  202137.0     1 MONTEVIDEO       UYMO
2     11928  722385.9     2    ARTIGAS       UYAR
3      4536  404534.8     3  CANELONES       UYCA
4      6106  431056.4     5    COLONIA       UYCO
5     11643  821325.4     6    DURAZNO       UYDU
6     10417  669409.9     8    FLORIDA       UYFD

• nombre.prj (opcional): sistema de referencia de coordenadas asociado a los datos.

crs_depart <- st_crs(depart_uy)
crs_depart$input

[1] "WGS 84 / UTM zone 21S"

• nombre.shx: almacena el índice de las entidades geométricas.

• nombre.sbn y .sbx (opcional): almacena índice espacial.

Como vimos en el ejemplo, read_sf() recibe la ruta a la carpeta o a alguno de los archivos de la capa de interés (en caso de tener varias capas en una carpeta podemos indicar la ruta hasta alguno de los archivos): read_sf("datos/ine_depto/ine_depto.shp") o utilizar el argumento layer: read_sf("datos/ine_depto", layer = "ine_depto")).

Geometrías

El segundo aspecto nuevo, es la columna geometry en los datos,

depart_uy$geometry

Geometry set for 20 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 366582.2 ymin: 6127919 xmax: 858252.1 ymax: 6671738
Projected CRS: WGS 84 / UTM zone 21S
First 5 geometries:

ésta puede almacenar distintos tipos de datos vectoriales.

Los tipos van aumentando en complejidad, el más simple serían puntos, cuya unión forma líneas, con estos dos formas polígonos. Algunas de las geometrías más comunes están representadas en la siguiente imagen:

Figure 1: Lovelace, Nowosad, y Muenchow (2019) (CC BY-NC-ND 4.0).

Tu turno

🗺 Grafica los polígonos de barrios de Montevideo (los archivos se encuentran en datos/ine_barrios).

🌎 Lee el archivo en datos/uptu_cultura y agrega las ubicaciones contenidas a la misma gráfica que los polígonos.

🌐 Describe la diferencia entre la geometría de los barrios de Montevideo y del archivo de cultura.

🤸 Un ejemplo más

Notemos que podemos utilizar las operaciones usuales para transformar tablas de datos, por ejemplo, si queremos crear una columna de tipo de establecimiento usamos la función mutate de dplyr.

glimpse(cultura_mv)

Rows: 398
Columns: 3
$ NOMBRE    <chr> "TEATRO ALFREDO MORENO", "TEATRO ALIANZA URUGUAY-E…
$ DIRECCION <chr> NA, "PARAGUAY 1217", "SAN JOSE 1426", "COLONIA 187…
$ geometry  <POINT [m]> POINT (582125 6138279), POINT (573785.3 6136…

cultura_mv <- cultura_mv %>% 
  mutate(tipo_establecimiento = case_when(
    str_detect(NOMBRE, "TEATRO") ~ "teatro", 
    str_detect(NOMBRE, "MUSEO") ~ "museo",
    str_detect(NOMBRE, "BIBLIOTECA") ~ "biblioteca",
    str_detect(NOMBRE, "CINE") ~ "cine",
    TRUE ~ "otro"
  ))

Usamos la función case_when() para categorizar los puntos de interés culturales de acuerdo al tipo, queremos esta variable para mostrar la información en un mapa.

glimpse(cultura_mv)

Rows: 398
Columns: 4
$ NOMBRE               <chr> "TEATRO ALFREDO MORENO", "TEATRO ALIANZ…
$ DIRECCION            <chr> NA, "PARAGUAY 1217", "SAN JOSE 1426", "…
$ geometry             <POINT [m]> POINT (582125 6138279), POINT (57…
$ tipo_establecimiento <chr> "teatro", "teatro", "teatro", "teatro",…

Veamos que tenemos disponibles las opciones típicas de ggplot como facet_wrap().

ggplot() +
  geom_sf(data = barrios_mv) +
  geom_sf(data = cultura_mv, color = "red", alpha = 0.8, size = 0.5) + 
  labs(subtitle = "Establecimientos culturales en Montevideo", 
       color = "tipo") +
  theme_void() + 
  facet_wrap(~tipo_establecimiento)

📃 Crear uno objeto sf de una tabla de datos

Imaginemos que tenemos una tabla de datos georreferenciados, pueden ser datos colectados con gps en una encuesta o ubicaciones obtenidas de google maps. Lo podemos almacenar en un data.frame que podemos a su vez convertir en un objeto de tipo sf, lo que debemos hacer es indicar cuáles son las columnas que incluyen información geográfica.

• 🖍 Utiliza google maps para ubicar un lugar especial para ti y agrégalo a este documento.

library(googlesheets4)
url <- ("https://docs.google.com/spreadsheets/d/1h7n3dpjqGofWYtKf5ADCDAr8nT5U5UalMFq9nS3G_ew/edit#gid=0")
gs4_deauth()
colab <- read_sheet(url)
glimpse(colab)

Rows: 3
Columns: 3
$ nombre <chr> "Teresa", "María", "casa"
$ long   <dbl> 18.93790, 19.36930, 18.95584
$ lat    <dbl> -99.23250, -99.12646, -99.22654

Por ahora colab es simplemente un data.frame,

class(colab)

[1] "tbl_df"     "tbl"        "data.frame"

Lo podemos convertir a un objeto de clase sf con la función st_as_sf, indicando cuales son las columnas con información geográfica. De manera opcional podemos indicar también la información de la proyección (Google Maps utiliza la 4326 para reportar).

colab_longlat <- colab %>% 
  st_as_sf(coords = c("lat", "long")) %>% 
  st_sf(crs = 4326)

class(colab_longlat)

[1] "sf"         "tbl_df"     "tbl"        "data.frame"

Y ahora lo podemos graficar con ggplot y geom_sf().

ggplot(colab_longlat) +
  geom_sf(aes(color = nombre))

📐 CRS y Proyecciones

Como vimos con nuestros primeros mapas, los objetos espaciales suelen tener asociado un valor en sistema de referencia de coordenadas (CRS), en esta sección explicamos a qué se refiere y porque es importante.

La Tierra tiene forma elipsoidal y para ubicar puntos en la Tierra debemos especificar un modelo de elipse, un punto de origen y dirección de los ejes, y una estrategia para proyectar a un plano.

1. Elipse: Describe la forma de la Tierra, es una aproximación que no ajusta a la perfección. Existen distintos elipsoides, algunos están diseñados para ajustar toda la Tierra (WGS84, GRS80) y algunos ajustan únicamente por regiones (NAD27). Los locales son más precisos para el área para la que fueron diseñados pero no funcionan en otras partes del mundo.

rgdal::projInfo(type = "ellps") %>% 
  rmarkdown::paged_table()

2. Datum: Define un punto de origen para el sistema de coordenadas de la Tierra y la dirección de los ejes. El datum define el elipsoide que se usará pero el sentido contrario no es cierto.

1. Las distancias deben tener unidades (metros en el ej.).
2. Se necesita un punto de referencia (datum).

Los conceptos anteriores son la base de un sistema de coordenadas, es importante recalcar que no todos los datos geoespaciales usan el mismo sistema de coordenadas. Y para ser representados en un mismo mapa debemos cuadrarlos.

3. Proyección: Nos permiten desplegar el elipsoide en un espacio de dos dimensiones, es decir pasar de longitud y latitud en el elipsoide a coordenadas de un mapa, esto facilita ciertas operaciones (por ejemplo. las longitudes y los ángulos son constantes en dos dimensiones, en contraste con una esfera donde las líneas de longitud convergen en los polos), sin embargo, conlleva necesariamente una distorsión. La estrategia usual para proyectar es utilizar superficies de desarrollo y después aplanar la superficie.

Figure 2: Comparación entre proyecciones tangente y secante cilíndirca, cónica y de azimut con licencia CC-BY 4.0

Es posible trabajar con datos no proyectados (longitud/latitud), a estos se les determina coordenadas geográficas.

Las distorsiones resultan en pérdidas de información que pueden ser en área, forma, distancia o dirección. Por ejemplo Mercator preserva dirección y es útil para navegación:

Show code

states <- map_data("state")
usmap <- ggplot(states, aes(x=long, y=lat, group=group)) +
  geom_polygon(fill="white", colour="black")
usmap + coord_map("mercator")

Por su parte Azimuthal Equal Area preserva area pero no dirección:

Show code

usmap + coord_map("azequalarea") 

En particular en investigación es común usar Universal Transverse Mercator (UTM) porque preserva ángulos y dirección (por lo que es común en sistemas de navegación). Para minimizar la distorsión se divide la Tierra en 60 regiones y utiliza una proyección (de secantes) en cada una. Esta es la proyección más común para Uruguay.

Figure 3: Husos y Zonas UTM con licencia CC-BY 3.0

🗃 Sistemas de referencia de coordenadas (CRS)

Como vimos arriba, los conjuntos de datos geoespaciales tienen asociado un sistema de referencia de coordenadas (y proyección si es el caso), Hay varias maneras de hacer referencia al sistema asociado a un conjunto de datos, una notación estándar que se utiliza para describir un CRS es proj4string y se ven como:

st_crs(barrios_mv)$proj4string # cartesianas

[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

st_crs(colab_longlat)$proj4string # geográficas

[1] "+proj=longlat +datum=WGS84 +no_defs"

También se puede usar el código ESPG asociado, los códigos EPSG. Los códigos EPSG (European Petroleum Survey Group) son códigos de 4 a 5 dígitos que representan definiciones de CRS:

st_crs(32721)$proj4string

[1] "+proj=utm +zone=21 +south +datum=WGS84 +units=m +no_defs"

Hay más maneras de documentar el sistema de coordenadas, sin embargo las dos mencionadas son muy comunes.

Tu turno

🗾 En un editor de texto abre el archivo .prj de alguno de los conjuntos de datos geoespaciales que hemos utilizado.

🗾 ¿Cuáles hemos usado hasta ahora?

Muchas veces es necesario transformar nuestros datos a otro sistema de coordenadas geográficas, por ejemplo, si queremos trabajar con conjuntos de datos con distinto CRS, o para utilizar ciertas herramientas como leaflet.

📌 Interactivos

Para hacer mapas interactivos usaremos el paquete leaflet (Cheng, Karambelkar, y Xie 2021):

library(leaflet)

Leaflet es una librería de JavaScript para realizar mapas interactivos, el paquete de R nos permite utilizarla desde R.

🖼 Puntos culturales con Leaflet

Leaflet recibe datos geoespaciales con CRS 4326 (al igual que las coordenadas que copiamos de GoogleMaps). Por ellos si queremos graficar los puntos culturales debemos comenzar reproyectando:

cultura_mv_longlat <- st_transform(cultura_mv, 4326)

Ahora graficamos, la sintaxis de leaflet es similar a la de ggplot, pero usamos %>% en lugar de +:

leaflet(data = cultura_mv_longlat) %>% 
  addTiles() %>% 
  addCircles()

La función addTiles() agrega las tejas del fondo. En este ejemplo agregagmos marcadores que muestran la ubicación de los sitios culturales de Montevideo, podemos añadir información de variables en los datos, por ejemplo, el nombre de la atracción. Y con addCircles() agregamos la geometría de puntos representados con círculos.

leaflet(data = cultura_mv_longlat) %>% 
  addTiles() %>% 
  addMarkers(popup = ~NOMBRE)

Con leaflet hay mucha flexibilidad, podemos meter html a las etiquetas o definir nuestros propios íconos, en el sitio se explica claramente con ejemplos.

rLadiesIcon <- makeIcon(
  iconUrl = "https://raw.githubusercontent.com/rladies/starter-kit/master/logo/R-LadiesGlobal_CMYK_offline_LogoOnly.svg",
  iconWidth = 38, iconHeight = 38
)

leaflet(data = colab_longlat) %>% 
  addTiles() %>% 
  addMarkers(popup = ~nombre, icon = rLadiesIcon) 

También podemos agregar objetos espaciales de otros tipos, para polígonos se utiliza addPolygons().

barrios_longlat <- barrios_mv %>% 
  st_transform(crs = 4326)
leaflet(data = barrios_longlat) %>% 
  addTiles() %>% 
  addPolygons() 

Tu turno

🌍 ¿Qué ocurre si utilizas leaflet con los datos no transformados, es decir si utilizas colab en lugar de colab_longlat?

🌍 Agrega al siguiente código la ubicación de los puntos de interés de Montevideo.

leaflet(data = colab_longlat) %>% 
  addTiles() %>% 
  addMarkers(popup = ~nombre, icon = rLadiesIcon) 

🌍 Experimenta con otros fondos, utiliza la función addProviderTiles() y elige un fondo de aquí.

leaflet(data = colab) %>% 
  addProviderTiles(providers$Stamen.Watercolor) %>% 
  addMarkers(popup = ~nombre, icon = rLadiesIcon) 

➖Extra: Operaciones

Las operaciones con datos geoespaciales van más allá del alcance de este taller, pero aquí hay algunos ejemplos de lo que se puede hacer.

• Área de polígonos (si no incluyera en los datos).

barrios_mv %>% 
  mutate(area = st_area(geometry))

Simple feature collection with 62 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 551982.7 ymin: 6133499 xmax: 589227.2 ymax: 6159810
Projected CRS: WGS 84 / UTM zone 21S
# A tibble: 62 x 5
   AREA_KM NOMBBARR   NROBARRIO                      geometry     area
 *   <dbl> <chr>          <int>            <MULTIPOLYGON [m]>    [m^2]
 1   2.11  Ciudad Vi…         1 (((573146.3 6137145, 573148.… 2297414…
 2   1.30  Centro             2 (((574497.2 6137037, 574506.… 1294508…
 3   0.694 Barrio Sur         3 (((574564 6136237, 574570.5 …  726727…
 4   2.28  Cordon             4 (((574397.4 6137451, 574382.… 2278114…
 5   0.798 Palermo            5 (((575149.1 6136881, 575144 …  832445…
 6   0.766 Parque Ro…         6 (((576530.7 6135354, 576532.…  778529…
 7   2.73  Punta Car…         7 (((576200.9 6135810, 576200.… 2818381…
 8   4.16  Buceo              9 (((580161.9 6137234, 580167.… 4164074…
 9   3.38  Pque. Bat…        10 (((578338.6 6137972, 578352.… 3371317…
10   3.50  Malvin            11 (((580642.4 6138757, 580590.… 3502102…
# … with 52 more rows

• En que barrio están los museos?

barrios_cultura_mv <- barrios_mv %>%
    st_join(cultura_mv, join = st_intersects) # usualmente st_within es más apropiado 

barrios_cultura_mv %>% 
  st_drop_geometry() %>% 
  group_by(NOMBBARR) %>% 
  summarise(n_lugares = sum(!is.na(NOMBRE))) %>% 
  arrange(-n_lugares) 

# A tibble: 62 x 2
   NOMBBARR                 n_lugares
   <chr>                        <int>
 1 Centro                          81
 2 Cordon                          57
 3 Ciudad Vieja                    50
 4 Pque. Batlle, V. Dolores        24
 5 Buceo                           14
 6 Palermo                         13
 7 Punta Carretas                  13
 8 Parque Rodo                     10
 9 Pocitos                         10
10 Aguada                           8
# … with 52 more rows

• Distancia entre puntos.

colab_longlat %>% 
  st_distance()

Units: [m]
          [,1]     [,2]      [,3]
[1,]     0.000 49039.48  2082.469
[2,] 49039.476     0.00 46963.494
[3,]  2082.469 46963.49     0.000

La cheatsheet del paquete sf describe operaciones geométricas disponibles, e incluye útiles esquemas.

Licencia

Texto y gráficas con licencia Creative Commons Attribution CC BY 4.0, a menos que se indique lo contrario o tengan otra fuente (indicado en la nota). Código en https://github.com/tereom/taller-mapas-mv.